Geometry, mechanics and electronics of singular structures and wrinkles in graphene 
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As the thinnest atomic membrane, graphene presents an opportunity to combine geometry, elas- 
ticity and electronics at the limits of their validity. The availability of reliable atomistic potentials 
for graphene allows unprecedented precise simulations of the mechanical response of atomic mem- 
branes. Here we describe the transport and electronic structure in the neighbourhood of conical 
singularities, the elementary excitations of the ubiquitous wrinkled and crumpled graphene that 
occur in non-epitaxial suspended samples where shear stresses are unavoidable. We use a combina- 
tion of atomistic mechanical simulations, analytical geometry and transport calculations in curved 
graphene, and exact diagonalization of the electronic spectrum to calculate the effects of geome- 
try on electronic structure, transport and mobility in suspended samples. We also point out how 
the geometry-generated pseudo-magnetic/electric fields might disrupt Landau quantization under a 
magnetic field. 



Graphene wrinkles easily and often [1]. This effect 
is most clearly observed in samples obtained from ex- 
foliation of graphite, and subsequent deposition onto a 
substrate [2], or in chemically derived oxides [3]. Since 
graphene is an atomically thin membrane, it is impossi- 
ble to lay a shear- free sheet of it onto a flat surface, as it 
sticks almost immediately to a substrate — such as the 
edges of a trench via van der Waals interaction — , and 
the substrate is itself rarely, if ever, flat [4-6], so that 
perfect shear-free conformations are not possible. In ad- 
dition, recently advanced techniques to grow graphene on 
metallic surfaces clearly show widespread wrinkling aris- 
ing from thermal expansion mismatch between graphene 
and the host substrate [7, 8]. These boundary deforma- 
tions acting on graphene lead to wrinkling because of the 
nearly negligible threshold for ID and 2D buckling in- 
stabilities in thin plates; the bending rigidity scales with 
the cube of the thickness so that a thin membrane cannot 
support even arbitrarily small shear or compression with- 
out wrinkling on scales large compared to its thickness 
[9]. 

However, for all its flexural limpness, graphene ex- 
hibits the largest in-plane Young's modulus [10] and, 
though easy to bend, is extremely hard to stretch. This 
geometry-induced separation of the energy scales for thin 
membranes implies that they try to respond to shear by 
bending isometrically almost everywhere [11]. However, 
except in very limited cases corresponding to developable 
deformations, bending alone cannot accommodate the 
state of stress or the boundary conditions imposed by 
the geometry. This conflict is resolved naturally by lo- 
cal membrane stretching by an amount sufficient to just 
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accommodate the imposed geometric and physical con- 
straints, so that regions of in-plane strain are restricted to 
vanishingly small areas distributed throughout the sys- 
tem. A simple example is seen in a thin sheet of paper 
which is very resistant to stretching (it actually tears be- 
fore we can stretch it), but bends easily; when a piece of 
crumpled paper is straightened out, we see flat areas con- 
nected by a network of ridges that meet at sharp vertices: 
the highly localized scars where the sheet is plastically 
deformed. The peaked structures constitute the basic 
element of the entire ridge network, and serve to focus 
large strains and energy densities. They are ubiquitous 
in thin films that are strongly deformed in such instances 
as drapes [11], skin winkles, etc., and are termed devel- 
opable cones or conical singularities (CSs); their outer 
form has been geometrically and mechanically character- 
ized in terms of a theory for the inextensional deforma- 
tions of thin sheets Ref. [12, 13]. In particular, there is 
a simple universal analytic expression for their geometry 
as a function of the boundary and/or stress conditions 
on the sheet far from the nearly singular tip where the 
effects of stretching are concentrated. 

Exfoliated graphene and suspended samples [2] natu- 
rally exhibit CSs as shown in Fig. 1(a), which are of par- 
ticular interest in the quest for ultimate electronic mo- 
bility [14, 15] and non-trivial interaction effects [16, 17]. 
Here we build on our knowledge of the structure and 
mechanics of CSs to study the influence of these ubiqui- 
tous objects on electronic transport in graphene. Un- 
like in most solid-state materials, flexural and planar 
deformations couple to electrons in graphene in a pe- 
culiar way. The honeycomb lattice implies that the ef- 
fective electronic excitations of the system are two di- 
mensional massless Dirac fermions [18]. The geometry of 
and strain in the lattice then couples to these excitations 
through both effective gauge fields, and local scattering 
potentials that follow the local curvature and thus affect 
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the electronic structure [19, 20], opening the prospect 
for strain-engineered electronic devices [21, 22]. CSs are 
also present in buckled nanotubes, where they have been 
shown to significantly alter transport characteristics [23] . 



I. CONICAL SINGULARITIES 

Graphene does behave like a thin plate under stress, 
even at the atomistic level; when sheared biaxially, 
and afterwards allowed to relax via Molecular Dynam- 
ics (MD) [24], Fig. 1(c) shows the relaxed configuration, 
which exhibits the classical Miura-ori like ridge pattern of 
2D buckling [25], with the CSs arising at the intersection 
of the ridges Fig. l(e,f). 

In cylindrical coordinates, the displacement field asso- 
ciated with the CSs, a generalized cone, reads 0) = 
r-ro = Up{p, e)up^ue{p, 0)ue^C{pi 0)z, where C(p, 0) = 
pip{0). The solution for ue{p^O)^Up{p^O)^ and il){0) is ob- 
tained by solving the equations of equilibrium for the 
finite bending of a plate with the constraint of inextensi- 
bility, i.e. that there is no in-plane strain (7^^ = 0). The 
vertical displacement is then given by (Refs. [12, 24]) 

C = m = £0(1^1 - ^1) + e^iemei - \e\) ; (1) 

where £ characterizes the angle of the enveloping cone, 
and both Oi ~ 70° and ^(^) are universal. This indepen- 
dence of the shape on any material parameters and scale, 
together with the Cauchy-Born hypothesis, allows us to 
describe conical singularities and wrinkling in graphene 
by applying the deformation field 6) to all atoms in 
the lattice. The resulting shape of the lattice is the one 
shown in Fig. 1(e), with the main effects arising from 
curvature. Since n(p, 0) is constructed so that there 
is no in-plane strain; however some localized stretch- 
ing strain is concentrated in the neighbourhood of the 
apex which will relax naturally in a MD simulation as a 
consequence of the large but finite stretching rigidity of 
graphene so that, even after relaxation, all inter- atomic 
distances are strongly peaked about the natural lattice 
spacing a = ao = 1.42 A, as can be seen in Fig. 1(d), with 
a spread of 2% for the values of e of interest here. This 
is just a refiection of the relative inextensibility of the 
in-plane a bonds, which leads to a blunting of the apex 
but is of little significance elsewhere. Since the relaxed 
structure shows strain > 1% for a dozen of atoms only, 
and very near the apex, we shall neglect it altogether. 



II. EFFECTIVE MODEL 

To understand how the electronic properties respond 
to this deformation, we note that the relevant physics 
occurs in the p^-derived tt bands of graphene; curvature 
causes re- hybridization of these orbitals [26], hindering 
or favouring wavefunction overlap, and thus perturbs the 
electronic kinetic energy. This affects both the tt band 



sub- system and hybridizes the Pz and the sp'^ sub-bands, 
which are otherwise orthogonal. As a first step we shall 
neglect this latter effect, which mostly shifts the chemical 
potential, and focus only on the tt bands. 

Within the tight-binding approximation, the band- 
structure is then determined by the effective Hamiltonian 

^ = E ^^^-4^. + E + H. c. , (2) 

where the two contributions come from first and second 
neighbours, and t['j = Vpp^ is the two centre Slater- 
Koster overlap integral [27], which has to be calculated 
now for all pairs of neighbours, taking into consideration 
the full geometry of the deformed lattice. To do this, we 
introduce the unit normal at every point of the surface, 
n(p, ^), so that for two atoms separated by an arbitrary 
distance d = Ri — Rj^ straightforward rotation of the Pz 
orbitals and Slater-Koster tables tell us that the overlap 
integral is [28] 

Uj = Vpp^ rii • rij + (Vppcj - Vppj,) {ui • d){nj ■ d) . (3) 

Since the surface is completely parametrized by the nor- 
mal displacement field, we may use the geometry of the 
developable cone to obtain the normals, distances and the 
hopping tij among any two atoms, noting that the under- 
lying metric remains Euclidean. To make progress ana- 
lytically we assume, with Harrison [29], that d?Vppx{d) = 
doVppx{do)^ so that on solving the Gauss equations we 
obtain tij{d) = t^j{do) + 5tij^ with: 

% « -V,%l\ido ■ V)VC|' + I [(do • V)2C]^ (4) 

Vi = Vpp^/3 — In the low energy approximation, 

we may then describe (2) by the effective Dirac Hamilto- 
nian: H ^ vfCt. [p - ^A{r)] + [3^0 + ^{r)] cr^ in each 
valley of the Brillouin zone. Then the effective gauge 
field w4(r) and the local potential ^{r) depend, respec- 
tively, on the perturbations of the nearest neighbour, and 
next-nearest neighbour hopping Ref. [18] via 

A,-iAy = J2SUr)e"'-^, $ = ^,54(r)e''=-^. (5) 

n A 

Substituting (4) in (5), one obtains 

^{r) = aTT^[d'djC] - (3 det[d'djC], (6) 

with a = 9alV^p^/8 + 27alVppa/32 ^ 1.5eVA^ and 

/3 = SalV^p^ + 9^^o^pV/8 - 3eVA^ [30]. We recah that 
d^djC, ~ K^j is the curvature tensor of our conical surface 
and, since H = and detK'^j are the local mean 

and Gaussian curvatures, it follows that ^ is entirely de- 
termined by the cone geometry. Moreover, since CSs are 
developable surfaces, the Gaussian curvature vanishes ev- 
erywhere, so that ^{r) = aao(V^C)^- The gauge field 
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A. is also given in terms of products of d^djC^ but we 
shall not write it explicitly since this potential couples 
to the electric current, and therefore does not contribute 
in leading order for scattering and transport when time- 
reversal symmetry is preserved. However $ leads to an 
electrostatic potential that is felt by the Dirac electrons 
and thus contributes directly to the resistivity. 

III. TRANSPORT 

We now consider the contribution of the CSs for the 
momentum relaxation time in the Boltzmann formal- 
ism. In the Born approximation, the scattering rate 
is given by S{k,k') = 27T/h\Vk,k'\^S{Ek - E^'), with 
Vfc,fc/ = ^k-k'[^ + exp(i^fc - i0fc/)]/2, and is the 
Fourier transform of the local potential (6): ^(r) = 
aa^[(lj{0) + {9)] / r'^ ^ which is of course directly related 
to the cone geometry. This potential is unusual for two 
reasons: it is anisotropic on account of (1) and decays 
in space as oc so that it is beyond the supercrit- 

ical threshold for Dirac fermions in 2D [31]. Were it 
not for the natural lattice regularization at r ~ 0, such 
potential would lead to an unbound spectrum of discrete 
states. This effect is also blunted the mechanical, stretch- 
induced relaxation observed near the apex in Fig. l(e,f). 
The result is a short range potential with a finite number 
of bound states (unlike the Coulomb case where the long 
range 1/r tail begets an infinite spectrum of resonances, 
even after regularization), so that CSs therefore scatter 
as short range, anisotropic potentials. 

The decay in the potential leads to an infrared 
divergence in with a leading order isotropic contribu- 
tion ^q ^ — lOae^ log(gro), all anisotropy being hidden 
in the subleading terms, with the regularization distance, 
ro of the order of the lattice spacing, refiecting the relax- 
ation in the neighbourhood of the apex. Then, the CSs 
scatter primarily as an isotropic potential, and the 
scattering time for the potential / can be calculated 
exactly, and reads 

7(M " V^ ' 0' 0' 0' -2 ; ' ^ ' 

where G is a Meijer function [32], the density of scat- 
terers, and relates to the carrier density via kp = irUe. 
Then the longitudinal conductivity follows from Eq. (7) 
and yields 

which is only relevant in the regime < kr^ < 1 shown 
in Fig. 2. We see that the conductivity is essentially 
linear in electron density throughout most of the region 
of interest, except for the logarithmic singularity around 
the Dirac point, where it grows sub-linearly. The cor- 
responding approximate mobility is /i ~ ^vf^^^o ^j^^i 



when is replaced by the corresponding parameter for 
CSs {uo 67ae^/7r eVA ) one obtains the mobility for a 
sea of uncorrelated CSs as /i ^ 10^^ro/(ni£:^j cm^/(Vs). 
Substituting the parameter values ro ~ 5 A and ~ 
lO^^cm-^ results in /i - lO^/e"^ cmV(Vs). The e"^ de- 
pendence refiects a strong sensitivity to the aperture of 
the enveloping cone of each CS, but given that £ < 0.5, it 
causes relatively small scattering. This effect should thus 
be more important in high- mobility suspended samples, 
where the CSs can become a limiting factor in carrier 
mobility. 



IV. ELECTRONIC SPECTRUM 

Although the gauge fields A. are not expected to con- 
tribute to transport at leading order, they do infiuence 
the electronic spectrum. In fact, since they arise from 
perturbations to nearest neighbour hopping, they might 
cause considerable fictitious magnetic fields [22]. To ad- 
dress this at the level of the lattice, we have calculated 
the electronic structure associated with the full tight- 
binding Hamiltonian (2) in the presence of a single un- 
relaxed CS. The local density of states [LDOS, Nr{E)] 
for representative parameters is shown in Fig. 3. We see 
that CSs scatter strongly enough to create bound elec- 
tronic states as shown in Fig. 3(a,c) by the sharp peaks 
for states beyond the band edge, decaying rapidly away 
from the apex. In addition, the LDOS is very struc- 
tured at other energies within the band, signalling the 
formation of resonant states. This is more clearly visible 
in Fig. 3(c) where the sampling points lie in the region 
of higher curvature. In this case the LDOS curves show 
even stronger perturbation around the Dirac point. The 
local bandwidth is decreased, and the leading slope of 
Nr{E) around E ^ fluctuates, indicating renormalized 
Fermi velocities in the neighbourhood of the apex. In 
panels Fig. 3(d-e) we plot the real-space distribution of 
the LDOS at representative energies around the Dirac 
point, showing that the charge density is mostly local- 
ized in the apex, albeit with a "leak" along two rays that 
are at an angle ~ 24° with the axis of symmetry of the 
CSs, coinciding with the two zero curvature generators in 
the entire conical surface (1), and show clearly the role of 
curvature-induced confinement [24]. We also always see 
signals of "magnetic" oscillations around the Dirac point, 
as shown in Fig. 3 (a) [33]. Even though these studies are 
carried out with zero magnetic field, the presence of these 
locally varying fictitious fields is expected to influence 
Landau level quantization under a real magnetic field. 



V. DISCUSSION 

We have shown that CSs have potential to markedly 
affect electronic properties and transport in wrinkled 
graphene. They contribute a quasi linear-in-density con- 
ductivity, and might even limit mobility in suspended 
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samples at low temperatures, even in the dilute limit. 
We did not consider the anisotropy in the transport cal- 
culations, but it is likely to play a role in situations like 
Fig. 1(b), where a strong alignment might lead to coher- 
ent scattering. In suspended samples, curvature-induced 
disordered flux might be dominant and thus explain why 
the quantum Hall effect in 4-terminal suspended samples 
is so elusive. Our calculations suggest that it is possi- 
ble to engineer the electronic and transport properties in 
graphene by inducing and controlling CSs on demand 
using substrate shear, or by exploring the anomalous 
thermal expansion coefficient of graphene, as initiated in 
Ref. [2] . Finally, such strong impact of singular deforma- 
tions on the electronic system can pave new avenues of 
interplay between structure and electronics. Graphene, 



as seen, has clear and unprecedented advantages, insofar 
as both its mechanical response and electronic structure 
are easily and accurately modelled down to the atomic 
level. 
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Figure 1: Wrinkles in graphene and the origin of conical singularities, a, Folded graphene sheet resembling the draping of a 
textile, and graphene suspended over a trench, b [34]. Some regions with visible conical singularities (CSs) are highlighted, c. Relaxed 
atomistic profile of a portion of graphene under biaxial shear, displaying typical buckling ridges, d. Logarithmic histogram of the 
inter-atomic distances in the relaxed configuration f for two values of e. e,f. Profile of the CSs studied here (e = 0.1). The atomic 
positions are shown as generated by applying u{p,0) to all atoms, e, and after relaxation by MD, f. 
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Figure 2: DC conductivity under the influence of CS's. DC conductivity (8) versus the adimensional electron density x = (kpro)'^. 
The dashed line shows the best linear fit in the entire domain {y ^ 12. 9x). Top inset shows a^^da/dx, which is how the electronic 
mobility (/x) is frequently extracted experimentally. Bottom inset amplifies the region kp ^ 0, dominated by a log singularity. 
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Figure 3: LDOS around the apex of a CS. LDOS close to the apex of a CS with e = 0.3, in the regions I, a, and II, c, and at 
specific lattice positions, as specified in panel b. The curves in a,c were vertically shifted for clarity. In d-f we show the LDOS in real 
space for selected energies: E/Vpp^ = -0.1,0.3,0.4. More detailed plots are available in [24]. 
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Supplementary Material 



Parametrization of conical singularities 

The fundamental object associated with a crumpled 
or wrinkled graphene sheet is the conical singularity [12], 
for which there is an analytic solution that minimizes the 
elastic energy subject to the constraint of inextensional 
deformations. In cylindrical polar coordinates (p, 6>), the 
displacement field relative to the fiat state is given by 
u{p,0) = r - To = Up{p,e)up + ue{p,6)ue + C{p,0)z, 
where the vertical displacement C{p^^) = Pi^{^) is a gen- 
eralized cone. The azimuthal component of the displace- 
ment field tl^{0) is given by 

^{0) = ee{\o\ - Oi) + £^(^)e(^i - 1^1) . (9) 

where O(-) is the usual Heaviside function, the aperture 
angle of the enveloping cone is tt — 2tan~^£:, and the 
universal geometric function ^{9) is given by: 

sin 6^1 cos a6> — a sin a^i cos ^ 
sin Oi cos aOi — a sin aOi cos 9i ' 

where a ^ 3.8 and Oi ^ 70° are universal numbers. 
The general shape of the deformed surface is shown in 
Fig. 4(a). The domain \0\ < Oi corresponds to the region 
where the surface coincides with the cylindrical envelope 
cone. 

The local azimuthal curvature of the cone is [il^"{0) + 
il){0)]/ p. In Fig. 4(b) we present a contour plot of 
e[ilj"{0) + V^(6>)]. For the region \0\ > 6>i, the shape is 
not that of a perfect cone; instead the local mean curva- 
ture increases attaining a maximum at Oq ~ 47°, before 
decreasing with an infiection point at O2 ^ 24° and be- 
comes negative for \0\ < 62^ with a maximum at ^ = 0. 

Details of the atomistic calculations 

We start with a graphene disk of radius 5nm which 
is geometrically transformed into a conical shape follow- 
ing the inextensional analytic model of [12]. This is done 
by fixing the coordinates of the atoms of 2 rows near the 
free edge, while the remaining carbon atoms of the cone 
are allowed to relax via a classical molecular dynamics 
simulation at zero temperature, so that the analytic so- 
lution which is singular at the tip relaxes by stretching 
locally. The C-C interactions are modeled with a second 



generation reactive empirical bond order (REBO) poten- 
tial [35]. 



Details of the LDOS calculations 

We use a tt tight-binding Hamiltonian (2), including 
overlaps up to next-to-nearest neighbors on a real Honey- 
comb lattice of dimensions 2000 x 1000 Oq. The lattice was 
deformed to conform to the exact analytical profile of the 
conical indentation discussed in the main text. The apex 
of the cone lies at the center of the lattice. The hopping 
integrals Uj and t were determined from the analyti- 
cal result (3). The local density of states (LDOS) was 
calculated by recursivelly solving for the local Green's 
function, until convergence is achieved. The results for 
the LDOS, Nr{E)^ are hence numerically exact, up to 
a broadening of AE/Vpp^ = 0.01, employed for faster 
convergence. 



Real space distribution of LDOS 

The LDOS, Nr{E), has been calculated exactly, as dis- 
cussed above, at all lattice sites in the vicinity of the 
apex, for readability, the real-space images show only 
a portion - less than 1% in area — of the total lattice. 
Since our interest rest on the local electronic structure 
near the apex, we considered only one conical singularity. 
The large lattice size ensures that spurious edge effects 
do not interfere with the structure around the apex. 

In Fig. 5 we present Nr{E) at all those lattice posi- 
tions, r, for a constant energy, E. The energy is mea- 
sured in units of Vpp^. At each lattice site a disk, whose 
diameter is proportional to the magnitude of the LDOS, 
is drawn. The panels in Fig. 5 pertain to an interval of 
energies between E/Vpp^ = — 1.7 and E/Vpp^ = 2.1, at 
increments of 0.2. This interval encompasses the two van 
Hove singularities, and the Dirac point which, for our 
choice of parameters, lies at E/Vpp^ = 0.3. 

For comparison, in Fig. 6 we present precisely the same 
analysis for a wider cone, with e = 0.1. In this case, 
there are no visible features, with the LDOS remaining 
greatly homogeneous, and following a variation with en- 
ergy much similar to the unperturbed lattice. Fig. 7 
shows an exaggerated case with e = 0.5. In this case 
relaxation effects should be stronger and the approxima- 
tion of neglecting a — n overlaps less warranted. Data is 
shown here for perspective only. 
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Figure 4: Geometry of a conical singularity, a, profile of a conical singularity with e — 0.2 seen from two different perspectives, 
b, Contour plot of £[(j:" {0) + revealing the curvature as a function of 0. The angles Oi ^ 70°, Oq ^ 46°, and O2 ~ 24° denote 

the generators along which the surface detaches from the enveloping cone, along which the curvature is maximal, and along which 
the curvature is null, respectively. 
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Figure 5: LDOS around the apex af a conical singularity. Each panel shows the distribution of the LDOS in real space around 
the apex of a cone with e — 0.3. Panels pertain to energies E /V^p^^ — —1.7, —1.5, —1.3, . . . , 1.9, 2.1. 
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Figure 6: LDOS around the apex af a conical singularity. Each panel shows the distribution of the LDOS in real space around 
the apex of a cone with e — 0.1. Panels pertain to energies E /V^p^^ — —1.7, —1.5, —1.3, . . . , 1.9, 2.1. 
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Figure 7: LDOS around the apex af a conical singularity. Each panel shows the distribution of the LDOS in real space around 
the apex of a cone with e — 0.5. Panels pertain to energies E /V^p^^ — —1.7, —1.5, —1.3, . . . , 1.9, 2.1. 



